! -----------------------------------------------
module okada92
! -----------------------------------------------
! calculates displacements due to a rectangular
! fault using Okada's formulation (1992).
! -----------------------------------------------
    implicit none

    ! structures 
    type disp
        real(8) :: x, y, z
    end type disp

    type fcomp
        real(8) :: u1, u2, u3
    end type fcomp

    type okadavars
        real(8) :: R          ! sqrt(xi^2+eta^2+p^2)
        real(8) :: yt, dt, ct ! y-tilde, d-tilde, c-tilde
        real(8) :: theta      ! atan(xi*eta/(q*R))
        real(8) :: X          ! xi^2+q^2
        real(8) :: lnRpXi     ! dlog(R+xi)
        real(8) :: lnRpEta    ! dlog(R+eta)
        real(8) :: X11, Y11, X32, Y32, Z32
    end type okadavars

    ! constants
    real(8), private, parameter :: pi = 3.1415926535897932d0

    ! functions
    public  :: okada92_disp
    private :: Chinnerys
    private :: fxyz2fcomp, fxyz2fcompC, add_fcomp
    private :: f_A_strike, f_B_strike, f_C_strike
    private :: f_A_dip, f_B_dip, f_C_dip
    private :: f_A_tensile, f_B_tensile, f_C_tensile
    private :: I1, I2, I3, I4, X11, Y11, X32, Y32, Z32, h
 
contains
    type(disp) function okada92_disp(x, y, z, a, c, dip, L, W, U1, U2, U3)
        real(8), intent(in) :: x, y, z  ! displacement at u(x,y,z)
        real(8), intent(in) :: a, c, dip, L, W, U1, U2, U3  ! fault geometry
        real(8) :: d, p, q  ! local coordinates on fault plane
        real(8) :: ux, uy, uz
        type(fcomp) :: ui_A, ui_Ah, ui_B, ui_C
        !real(8), parameter :: pi = 3.141592653589793238d0 -> moved to global in the module

        d = c - z
        p = y*dcos(dip) + d*dsin(dip)
        q = y*dsin(dip) - d*dcos(dip)

        ! displacement due to strike slip
        ui_A  = Chinnerys(f_A_strike, a, x, z, dip, p, q, L, W)
        ui_Ah = Chinnerys(f_A_strike, a, x,-z, dip, p, q, L, W)
        ui_B  = Chinnerys(f_B_strike, a, x, z, dip, p, q, L, W)
        ui_C  = Chinnerys(f_C_strike, a, x, z, dip, p, q, L, W)
        ux = U1/(2.d0*pi)*( ui_A%u1 - ui_Ah%u1 + ui_B%u1 + z*ui_C%u1)
        uy = U1/(2.d0*pi)*((ui_A%u2 - ui_Ah%u2 + ui_B%u2 + z*ui_C%u2)*dcos(dip) &
                        -  (ui_A%u3 - ui_Ah%u3 + ui_B%u3 + z*ui_C%u3)*dsin(dip))
        uz = U1/(2.d0*pi)*((ui_A%u2 - ui_Ah%u2 + ui_B%u2 - z*ui_C%u2)*dsin(dip) &
                        +  (ui_A%u3 - ui_Ah%u3 + ui_B%u3 - z*ui_C%u3)*dcos(dip))
        
        ! displacement due to dip slip
        ui_A  = Chinnerys(f_A_dip, a, x, z, dip, p, q, L, W)
        ui_Ah = Chinnerys(f_A_dip, a, x,-z, dip, p, q, L, W)
        ui_B  = Chinnerys(f_B_dip, a, x, z, dip, p, q, L, W)
        ui_C  = Chinnerys(f_C_dip, a, x, z, dip, p, q, L, W)
        ux = ux + U2/(2.d0*pi)*( ui_A%u1 - ui_Ah%u1 + ui_B%u1 + z*ui_C%u1)
        uy = uy + U2/(2.d0*pi)*((ui_A%u2 - ui_Ah%u2 + ui_B%u2 + z*ui_C%u2)*dcos(dip) &
                        -  (ui_A%u3 - ui_Ah%u3 + ui_B%u3 + z*ui_C%u3)*dsin(dip))
        uz = uz + U2/(2.d0*pi)*((ui_A%u2 - ui_Ah%u2 + ui_B%u2 - z*ui_C%u2)*dsin(dip) &
                        +  (ui_A%u3 - ui_Ah%u3 + ui_B%u3 - z*ui_C%u3)*dcos(dip))

        ! displacement due to strike slip
        ui_A  = Chinnerys(f_A_tensile, a, x, z, dip, p, q, L, W)
        ui_Ah = Chinnerys(f_A_tensile, a, x,-z, dip, p, q, L, W)
        ui_B  = Chinnerys(f_B_tensile, a, x, z, dip, p, q, L, W)
        ui_C  = Chinnerys(f_C_tensile, a, x, z, dip, p, q, L, W)
        ux = ux + U3/(2.d0*pi)*( ui_A%u1 - ui_Ah%u1 + ui_B%u1 + z*ui_C%u1)
        uy = uy + U3/(2.d0*pi)*((ui_A%u2 - ui_Ah%u2 + ui_B%u2 + z*ui_C%u2)*dcos(dip) &
                        -  (ui_A%u3 - ui_Ah%u3 + ui_B%u3 + z*ui_C%u3)*dsin(dip))
        uz = uz + U3/(2.d0*pi)*((ui_A%u2 - ui_Ah%u2 + ui_B%u2 - z*ui_C%u2)*dsin(dip) &
                        +  (ui_A%u3 - ui_Ah%u3 + ui_B%u3 - z*ui_C%u3)*dcos(dip))

        okada92_disp%x = ux
        okada92_disp%y = uy
        okada92_disp%z = uz
        return
    end function

    !real(8) function Chinnerys(f, a, x, z, dip, p, q, L, W)
    type(fcomp) function Chinnerys(f, a, x, z, dip, p, q, L, W)
        type(fcomp) :: f  ! function
        type(fcomp) :: f_A, f_Ah, f_B, f_C
        real(8) :: a, x, z, dip, p, q, L, W
        real(8) :: xi(4), eta(4), eps
        type(okadavars) :: v(4)
        integer :: i

        xi(1)  = x; xi(2)  = x;   xi(3)  = x-L; xi(4)  = x-L;
        eta(1) = p; eta(2) = p-W; eta(3) = p;   eta(4) = p-W;

        eps = epsilon(dip) ! epsilon of real(8).
        do i = 1, 4
            v(i)%R  = dsqrt(xi(i)*xi(i) + eta(i)*eta(i) + q*q)
            v(i)%yt = eta(i)*dcos(dip) + q*dsin(dip)
            v(i)%dt = eta(i)*dsin(dip) - q*dcos(dip)
            v(i)%ct = v(i)%dt + z
            v(i)%X  = dsqrt(xi(i)*xi(i) + q*q)
            v(i)%X11= X11(xi(i), v(i)%R)
            v(i)%Y11= Y11(eta(i),v(i)%R)
            v(i)%X32= X32(xi(i), v(i)%R)
            v(i)%Y32= Y32(eta(i),v(i)%R)
            v(i)%Z32= Z32(eta(i), v(i)%R, dip, q, z)
            ! spacial cases
            if (dabs(q) > eps) then
                v(i)%theta = datan2(xi(i)*eta(i), q*v(i)%R)
            else
                v(i)%theta = 0.d0 
            endif

            if (dabs(v(i)%R+xi(i)) > eps) then
                v(i)%lnRpXi = dlog(v(i)%R+xi(i))
            else
                v(i)%lnRpXi = -dlog(v(i)%R-xi(i))
            endif

            if (dabs(v(i)%R+eta(i)) > eps) then
                v(i)%lnRpEta = dlog(v(i)%R+eta(i))
            else
                v(i)%lnRpEta = -dlog(v(i)%R-eta(i))
            endif
        enddo

        !Chinnerys = f(a, xi(1), eta(1), z, dip, q, v(1))  &
        !          - f(a, xi(2), eta(2), z, dip, q, v(2))  &
        !          - f(a, xi(3), eta(3), z, dip, q, v(3))  &
        !          + f(a, Xi(4), eta(4), z, dip, q, v(4))

        f_A  = f(a, xi(1), eta(1), z, dip, q, v(1)) 
        f_Ah = f(a, xi(2), eta(2), z, dip, q, v(2))
        f_B  = f(a, xi(3), eta(3), z, dip, q, v(3))
        f_C  = f(a, Xi(4), eta(4), z, dip, q, v(4))
        Chinnerys%u1 = f_A%u1 - f_Ah%u1 - f_B%u1 + f_C%u1
        Chinnerys%u2 = f_A%u2 - f_Ah%u2 - f_B%u2 + f_C%u2
        Chinnerys%u3 = f_A%u3 - f_Ah%u3 - f_B%u3 + f_C%u3
        return
    end function

    ! for Part A and B
    type(fcomp) function fxyz2fcomp(fx, fy, fz, dip)
        real(8), intent(in) :: fx, fy, fz, dip
        fxyz2fcomp%u1 = fx
        fxyz2fcomp%u2 = fy*dcos(dip) + fz*dsin(dip)
        fxyz2fcomp%u3 =-fy*dsin(dip) + fz*dcos(dip)
        return
    end function

    ! for Part C
    type(fcomp) function fxyz2fcompC(fx, fy, fz, dip)
        real(8), intent(in) :: fx, fy, fz, dip
        fxyz2fcompC%u1 = fx
        fxyz2fcompC%u2 = fy*dcos(dip) - fz*dsin(dip)
        fxyz2fcompC%u3 =-fy*dsin(dip) - fz*dcos(dip)
        return
    end function

    type(fcomp) function add_fcomp(f1, f2)
        type(fcomp), intent(in) :: f1, f2
        add_fcomp%u1 = f1%u1 + f2%u1
        add_fcomp%u2 = f1%u2 + f2%u2
        add_fcomp%u3 = f1%u3 + f2%u3
        return
    end function

    type(fcomp) function f_A_strike(a, xi, eta, z, dip, q, v)
        real(8), intent(in) :: a, xi, eta, z, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: fx, fy, fz
        real(8) :: theta, R, lnRpEta, Y11

        theta = v%theta; R = v%R; lnRpEta = v%lnRpEta; Y11 = v%Y11;
        fx = 0.5d0*theta + 0.5d0*a*xi*eta*Y11
        fy = 0.5d0*a * q/R
        fz = (1.d0-a)/2.d0*lnRpEta - (a/2.d0)*(q**2.d0)*Y11
        f_A_strike = fxyz2fcomp(fx, fy, fz, dip)
        return
    end function

    type(fcomp) function f_B_strike(a, xi, eta, z, dip, q, v)
        real(8), intent(in) :: a, xi, eta, z, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: fx, fy, fz
        real(8) :: theta, R, yt, dt, Y11

        theta = v%theta; R = v%R; yt = v%yt; dt = v%dt; Y11 = v%Y11;
        fx = -xi*q*Y11 - theta - (1.d0-a)/a*I1(xi, eta, dip, q, v)*dsin(dip)
        fy = -q/R              + (1.d0-a)/a*yt/(R+dt)*dsin(dip)
        fz =  q**2.d0*Y11      - (1.d0-a)/a*I2(xi, eta, dip, q, v)*dsin(dip)
        f_B_strike = fxyz2fcomp(fx, fy, fz, dip)
        return
    end function

    type(fcomp) function f_C_strike(a, xi, eta, z, dip, q, v)
        real(8), intent(in) :: a, xi, eta, z, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: fx, fy, fz
        real(8) :: R, ct, Y11, Z32

        R = v%theta; ct = v%ct; Y11 = v%Y11; Z32 = v%Z32;
        fx = (1.d0-a)*xi*Y11*dcos(dip)                   - a*xi*q*Z32
        fy = (1.d0-a)*(dcos(dip)/R+2.d0*q*Y11*dsin(dip)) - a*ct*q/(R**3.d0)
        fz = (1.d0-a)*q*Y11*dcos(dip)  &
                      - a*(ct*eta/(R**3.d0)-z*Y11+(xi**2.d0)*Z32)
        f_C_strike = fxyz2fcompC(fx, fy, fz, dip)
        return
    end function

    type(fcomp) function f_A_dip(a, xi, eta, z, dip, q, v)
        real(8), intent(in) :: a, xi, eta, z, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: fx, fy, fz
        real(8) :: theta, R, lnRpXi, X11

        theta = v%theta; R = v%R; lnRpXi = v%lnRpXi; X11 = v%X11;
        fx = a/2.d0*q/R
        fy = theta/2.d0        + a/2.d0*eta*q*X11
        fz = (1-a)/2.d0*lnRpXi - a/2.d0*(q**2.d0)*X11
        f_A_dip = fxyz2fcomp(fx, fy, fz, dip)
        return
    end function

    type(fcomp) function f_B_dip(a, xi, eta, z, dip, q, v)
        real(8), intent(in) :: a, xi, eta, z, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: fx, fy, fz
        real(8) :: theta, R, dt, X11

        theta = v%theta; R = v%R; dt = v%dt; X11 = v%X11;
        fx = -q/R               + (1.d0-a)/a*I3(xi, eta, dip, q, v)*dsin(dip)*dcos(dip)
        fy = -eta*q*X11 - theta - (1.d0-a)/a*xi/(R+dt)*dsin(dip)*dcos(dip)
        fz =  q*q*X11           + (1.d0-a)/a*I4(xi, eta, dip, q, v)*dsin(dip)*dcos(dip)
        f_B_dip = fxyz2fcomp(fx, fy, fz, dip)
        return
    end function

    type(fcomp) function f_C_dip(a, xi, eta, z, dip, q, v)
        real(8), intent(in) :: a, xi, eta, z, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: fx, fy, fz
        real(8) :: R, yt, ct, dt, X11, Y11, X32

        R = v%R; yt = v%yt; ct = v%ct; dt = v%dt; X11 = v%X11; Y11 = v%Y11; X32 = v%X32;
        fx = (1.d0-a)*dcos(dip)/R - q*Y11*dsin(dip)  - a*ct*q/(R**3.d0)
        fy = (1.d0-a)*yt*X11                         - a*ct*eta*q*X32
        fz = -dt*X11              - xi*Y11*dsin(dip) - a*ct*(X11-q*q*X32)
        f_C_dip = fxyz2fcompC(fx, fy, fz, dip)
        return
    end function

    type(fcomp) function f_A_tensile(a, xi, eta, z, dip, q, v)
        real(8), intent(in) :: a, xi, eta, z, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: fx, fy, fz
        real(8) :: theta, R, lnRpXi, lnRpEta, X11, Y11

        theta = v%theta; R = v%R; lnRpXi = v%lnRpXi; lnRpEta = v%lnRpEta; X11 = v%X11; Y11 = v%Y11;
        fx = -(1.d0-a)/2.d0*lnRpEta - a/2.d0*q*q*Y11
        fy = -(1.d0-a)/2.d0*lnRpXi  - a/2.d0*q*q*X11
        fz =  theta/2.d0            - a/2.d0*q*(eta*X11+xi*Y11)
        f_A_tensile = fxyz2fcomp(fx, fy, fz, dip)
        return
    end function

    type(fcomp) function f_B_tensile(a, xi, eta, z, dip, q, v)
        real(8), intent(in) :: a, xi, eta, z, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: fx, fy, fz
        real(8) :: theta, R, dt, lnRpXi, lnRpEta, X11, Y11

        theta = v%theta; R = v%R; dt = v%dt; X11 = v%X11; Y11 = v%Y11;
        fx = q*q*Y11 - (1.d0-a)/a*I3(xi, eta, dip, q, v)*(dsin(dip))**2.d0
        fy = q*q*X11 + (1.d0-a)/a*xi/(R+dt)*(dsin(dip))*2.d0
        fz = q*(eta*X11+xi*Y11) - theta &
                                 - (1.d0-a)/a*I4(xi, eta, dip, q, v)*(dsin(dip))**2.d0
        f_B_tensile = fxyz2fcomp(fx, fy, fz, dip)
        return
    end function

    type(fcomp) function f_C_tensile(a, xi, eta, z, dip, q, v)
        real(8), intent(in) :: a, xi, eta, z, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: fx, fy, fz
        real(8) :: R, dt, yt, ct, X11, Y11, X32, Z32

        R = v%R; dt = v%dt; yt = v%yt; ct = v%ct; X11 = v%X11; Y11 = v%Y11; X32 = v%X32; Z32 = v%Z32;
        fx = -(1.d0-a)*(dsin(dip)/R+q*Y11*dcos(dip))  - a*(z*Y11-q*q*Z32)
        fy =  (1.d0-a)*2.d0*xi*Y11*dsin(dip) + dt*X11 - a*ct*(X11-q*q*X32)
        fz =  (1.d0-a)*(yt*X11+xi*Y11*dcos(dip))      + a*q*(ct*eta*X32+xi*Z32)
        f_C_tensile = fxyz2fcompC(fx, fy, fz, dip)
        return
    end function

    ! ---------- specific terms ----------!
    real(8) function I1(xi, eta, dip, q, v)
        real(8), intent(in) :: xi, eta, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: R, dt

        R = v%R; dt = v%dt;
        I1 = -xi/(R+dt)*dcos(dip) - I4(xi, eta, dip, q, v)*dsin(dip)
        return
    end function

    real(8) function I2(xi, eta, dip, q, v)
        real(8), intent(in) :: xi, eta, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: R, dt

        R = v%R; dt = v%dt;
        I2 = dlog(R+dt) + I3(xi, eta, dip, q, v)*dsin(dip)
        return
    end function

    real(8) function I3(xi, eta, dip, q, v)
        real(8), intent(in) :: xi, eta, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: R, dt, yt, lnRpEta
        
        R = v%R; dt = v%dt; yt = v%yt; lnRpEta = v%lnRpEta
        if (dabs(dcos(dip)) > epsilon(dip)) then
            I3 = yt/(dcos(dip)*(R+dt)) &
               !- (1.d0/(dcos(dip))**(2.d0))*(dlog(R+eta) - dsin(dip)*dlog(R+dt))
               - (1.d0/(dcos(dip))**2.d0)*(lnRpEta - dsin(dip)*dlog(R+dt))
        else
            !I3 = 0.5d0*(eta/(R+dt) + yt*q/(R+dt)**(2.0) -dlog(R+eta))
            I3 = 0.5d0*(eta/(R+dt) + yt*q/(R+dt)**2.d0 -lnRpEta)
        endif
        return
    end function

    real(8) function I4(xi, eta, dip, q, v)
        real(8), intent(in) :: xi, eta, dip, q
        type(okadavars), intent(in) :: v
        real(8) :: R, X, dt, yt

        R = v%R; X = v%X; dt = v%dt; yt = v%yt;
        if (dabs(xi) <= epsilon(xi)) then
            I4 = 0.d0
        elseif (dabs(dcos(dip)) > epsilon(dip)) then
            I4 = (dsin(dip)/dcos(dip)*xi/(R+dt)) &
               + (2.d0/(dcos(dip))**2.d0)*datan2(eta*(X+q*dcos(dip)) + X*(R+X)*dsin(dip), xi*(R+X)*dcos(dip))
        else
            I4 = 0.5d0*(xi*yt/(R+dt)**2.d0)
        endif
        return
    end function

    real(8) function X11(xi, R)
        real(8), intent(in) :: xi, R

        if (dabs(R+xi) > epsilon(R)) then
            X11 = 1.d0 / (R*(R+xi))
        else
            X11 = 0.d0
        endif
        return
    end function

    real(8) function Y11(eta, R)
        real(8), intent(in) :: eta, R

        if (dabs(R+eta) > epsilon(R)) then
            Y11 = 1.d0 / (R*(R+eta))
        else
            Y11 = 0.d0
        endif
        return
    end function

    real(8) function X32(xi, R)
        real(8), intent(in) :: xi, R

        if (dabs(R+xi) > epsilon(R)) then
            X32 = (2.d0*R + xi) / (R**3.d0*(R+xi)**2.d0)
        else
            X32 = 0.d0
        endif
        return
    end function

    real(8) function Y32(eta, R)
        real(8), intent(in) :: eta, R

        if (dabs(R+eta) > epsilon(R)) then
            Y32 = (2.d0*R + eta) / (R**3.d0*(R+eta)**2.d0)
        else
            Y32 = 0.d0
        endif
        return
    end function

    real(8) function Z32(eta, R, dip, q, z)
        real(8), intent(in) :: eta, R, dip, q, z
        Z32 = dsin(dip) / (R**3.d0) - h(dip, q, z)*Y32(eta, R)
    end function

    real(8) function h(dip, q, z)
        real(8), intent(in) :: dip, q, z
        h = q*dcos(dip) - z
    end function
end module
